#include <iostream>
#include <fstream>
#include <cmath>
#include <cstdlib>
#include <ctime>


using namespace std;

class SingleParticleSystem{

    public:
        SingleParticleSystem(double _temperature, long _times, double _delta){

            srand(time(NULL));

            temperature = _temperature;
            beta = 1.0 / temperature;

            v = 0;
            energy = 0.5 * v * v;
            energyAccumulation = 0;


            accept = 0;
            times = _times;

            delta = _delta;

        }



        int Run(int _times){
            for (long i = 0; i < _times; i++){
                step();
            }
            return 1;
        }

        void step(){
            double vTrial = v + delta * (2.0 * ran1() - 1.0);
            double energyTrial = getEnergyByVelocity(vTrial);
            double de = energyTrial - energy;
            if (de < 0 || exp(-1.0 * beta * de) > ran1()){
                accept++;
                v = vTrial;
                energy = energyTrial;
            }
            energyAccumulation += energy;
        }

        void getToEquilibrium(){
            Run((long)(times * 0.2));
            reset();
        }


        void reset(){
            energyAccumulation = 0;
            accept = 0;
        }

        ~SingleParticleSystem(){

        }


        double  getEnergyByVelocity(double _v){
            return 0.5 * _v * _v;
        }

        double ran1(){
            return 1.0 * rand() / RAND_MAX;
        }

        long getTimes(){
            return times;
        }


        double getEnergy(){
            return 1.0 * energyAccumulation / times;
        }

        double getVelocity(){
            return v;
        }

        double getAccept(){
            return 1.0 * accept / times;
        }


    private:
        double temperature;
        double beta;

        double v;
        double energy;
        double energyAccumulation;

        long accept;
        long times;
        double delta;



};


int main(){

    ofstream output("single_particle.txt");
    for (double t = 1.0; t < 300; t += 1.0){
        SingleParticleSystem s(t, 1000000, 4.0);
        s.getToEquilibrium();
        s.Run(s.getTimes());
        output << t << "\t" << s.getEnergy() << "\t" << s.getVelocity() << "\t" << s.getAccept() << endl;
    }
    output.flush();
    output.close();
    return 0;

}
